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1.0  ACTIVITY  SUMMARY 


This  is  the  Phase  I  final  report,  which  summarizes  our  following  R  &  D  activities: 

1.1  Development  of  a  BGK-NS  Solver  for  Magnetogasdynamics  in  Thermal  and  Chemical 
Equilibrium 

Because  of  the  time  limitations  of  Phase  I  and  the  large  amount  of  effort  invested  in  the  development  of 
the  gas-kinetic  flow  solver  and  the  adjoint  optimization  method  for  weakly  ionized  hypersonic  flows,  we 
have  focused  on  the  development  of  a  BGK-NS  solver  for  magnetogasdynamics  in  thermal  and  chemical 
equilibrium  in  Phase  I.  This  is  a  first  step  that  will  be  further  extended  in  Phase  II.  Section  2  discusses  the 
governing  equations  of  magnetogasdynamics  in  thermal  and  chemical  equilibrium  and  the  corresponding 
gas-kinetic  BGK-NS  solver. 

1.2  Development  of  an  Adjoint  Optimization  Method  for  Magnetogasdynamics  in  Thermal  and 
Chemical  Equilibrium 

Based  on  the  above  BGK-NS  solver,  we  have  developed  the  corresponding  discrete  adjoint  equations  for 
the  viscous,  low  magnetic  Reynolds  number  approximation.  These  equations  have  been  implemented 
into  a  preliminary  adjoint  solver  which  is  capable  of  producing  gradient  information  for  arbitrary  cost 
functions  (unlike  the  continuous  version  of  the  adjoint)  and  for  arbitrarily  large  numbers  of  design 
parameters  with  minimum  cost  (a  flow  solution  and  an  adjoint  solution  only.)  Details  of  this  work  are 
presented  in  Section  3. 

1.3  Hypersonic  Flow  Past  a  Cylinder:  Test  Case 

To  validate  the  gas-kinetic  BGK-NS  solver  and  the  adjoint  optimization  method,  the  hypersonic  flow  past 
cylinder  case  were  numerically  investigated  and  the  results  are  presented  in  Section  4. 

1.4  Development  of  BGK-Burnett  Solver  for  Gas  Dynamics  in  Thermal  and  Chemical 
Equilibrium 

In  addition,  we  will  develop  BGK-Bumett  solver  in  Section  5  to  explore  the  possibility  of  extending  the 
approach  beyond  the  continuum  regime.  In  order  to  validate  the  resulting  BGK-Bumett  scheme,  the  plane 
Poiseuille  flow  will  be  numerically  investigated. 

1.5  Development  of  Gas-Kinetic  Navier-Stokes  Solver  for  Weakly  Ionized  Hypersonic  Flows 

In  addition,  we  have  also  explored  the  possibility  of  extending  the  gas-kinetic  CFD  approach  to  weakly 
ionized  hypersonic  flows.  Following  an  11-species  air  model  and  the  two-temperature  model  used  in 
LAURA  [1-2],  a  preliminary  gas-kinetic  equivalence  of  LAURA  is  presented  in  Section  6. 
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2.0  BGK-NS  SCHEME  FOR  MAGNETOGASDYNAMICS 
IN  THERMAL  AND  CHEMICAL  EQUILIBRIUM 

Our  ultimate  goal  for  this  STTR  project  is  to  develop  an  adjoint  optimization  method  for  electromagnetic 
control  of  weakly  ionized  hypersonic  flows.  Due  to  the  time  limitations  of  Phase  I  and  the  large  amount 
of  initial  effort  involved  in  the  development  of  an  adjoint  optimization  method  based  on  a  gas-kinetic 
flow  solver,  however,  we  have  first  developed  a  gas-kinetic  flow  solver  and  an  adjoint  optimization 
method  for  magnetogasdynamics  in  thermal  and  chemical  equilibrium.  This  will  serve  as  the  foundation 
for  a  further  and  more  detailed  generalization  of  the  approach  to  weakly  ionized  hypersonic  flows  in 
Phase  n. 


2.1  Governing  Equations 

The  equations  governing  a  conducting  flow  in  a  magnetic  field  can  be  obtained  by  simply  coupling  the 
pre-Maxwell  equation  to  the  fluid  conservation  equations  through  the  momentum  and  energy  equations. 
The  fluid  conservation  equations  up  to  the  Navier-Stokes  order  of  the  Chapman-Enskog  expansion  are 


I-  W ) + St  W 'u‘  *  s‘ ‘p - T" )  -  0 

ot  dxJ 

—  {pa)+-^—[{ps  +  p)Uj  -U't‘j  -k^—]  =  0 
dt dx}  H  y  dxJ 


(2-1) 


Here  p ,  p ,  T ,  and  a  are  the  density,  pressure,  temperature,  and  the  total  energy  respectively.  U‘  is  the 
velocity  component  and  zij  is  the  shear  stress,  k  is  the  thermal  conductivity  coefficient. 


In  the  presence  of  a  magnetic  field,  a  conducting  fluid  can  generate  a  conduction  current  j‘ ,  which  may 
be  computed  by  the  pre-Maxwell  equation 


d  ,B 


r( — ) 


dxJ  p, 


(2-2) 


where  Bk  is  the  magnetic  field  and  pm  is  the  magnetic  permeability.  To  include  the  electromagnetic 
effects  into  (2-1),  one  should  add  the  Lorentz  force  F'em  =  sijk  jJBk  to  the  momentum  equation  and  Joule 

heating  Pem  =  E‘  j‘  to  the  total  energy  equation.  Here  E‘  is  the  electric  field,  which  may  be  computed  by 
combining  the  pre-Maxwell  equation  of  (2-2)  and  the  Ohm’s  law: 


E‘  =— —  aiJkUj  Bk 
cr 

where  cr  is  the  electrical  conductivity.  As  a  result,  equations  of  (2-1)  are  generalized  to 


(2-3) 
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(2-4) 


ip+ij{puJ)=0 

Jt  (pU‘ )+JL(pUiUJ+  8ij p -T'J)  =  sijk  jJ Bk 

| .(P£)+^T[(ps  +  P)UJ -U1^ -K^r)  =  ^—eiJkfUjBk 
dt  8xJ  dx1  cr 

A  conducting  fluid  in  an  imposed  magnetic  field  can  self  induce  a  magnetic  field.  If  the  influence  of  the 
self-induced  magnetic  field  is  not  negligible,  one  has  to  further  add  the  magnetic  induction  equation, 
which  is  derived  from  Faraday’s  law,  to  (2-4) 

—  p+—{pUj)  =  0 
dt H  dxJ 

—  (pu‘ )  +  —  (pUiUj  +Sijp-  Tij  )  =  eiJkjJBk 

dt  dx1  i  .  (2-5) 

^-(pe)  +  JL[(p£  +  p)Uj -U‘riJ  - k - siJk  j'lJJ Bk 
dt  dx1  dx1  cr 

|-5'-(l-^) -A-  (C/'^-iy^)  =  -(l-^)_^-{i[A(^i)- JL  (11)]} 

dt  dx1  dx1  cr  dx  pm  8xJ  pm 

For  multidimensional  cases,  there  is  another  Maxwell’s  equation  which  acts  as  a  constraint  to  the 
numerical  algorithm 


dB‘ 


dx‘ 


-  =  0 


(2-6) 


2.2  Gas-kinetic  BGK-NS  Scheme 

There  are  several  numerical  algorithms  available  in  the  literature  (e.g.,  [3-5])  for  the  solution  of  the 
governing  equations  of  (2-5)  based  on  the  continuum  CFD  approach.  In  this  project,  we  will  develop  a 
gas-kinetic  CFD  algorithm  for  the  solution  of  the  governing  equations  of  (2-5).  The  gas-kinetic  CFD 
approach  has  several  advantages  over  the  more  widely  used  continuum  CFD  approach: 

•  Similar  to  DSMC,  no  macroscopic  governing  equation  is  needed  before  the  construction  of  a  gas- 
kinetic  scheme.  Since  the  macroscopic  governing  equations  are  derived  from  the  gas-kinetic  theory, 
a  gas-kinetic  CFD  method  can  always  satisfy  the  corresponding  macroscopic  governing  equations 
automatically 

-  The  entropy  condition  is  automatically  satisfied 

-  The  scheme  is  positivity-preserving 

-  The  scheme  is  free  of  any  sonic  point  glitch 

-  The  scheme  is  free  of  odd-even  decoupling 

•  A  gas-kinetic  scheme  treats  the  inviscid  and  viscous  fluxes  as  a  single  entity.  This  makes  its 
extension  beyond  the  continuum  regime  much  easier.  Its  integration  with  DSMC  or  Direct 
Boltzmann  solvers  in  the  rarefied  flow  regime  is  also  much  more  straightforward  than  the 
continuum  CFD  approach. 

•  The  numerical  fluxes  given  by  a  gas-kinetic  scheme  are  for  all  speeds.  No  preconditioning  is 
needed  for  incompressible  or  hypersonic  flows. 

However,  theoretically  it  is  very  difficult  to  construct  an  equilibrium  state  and  a  single  kinetic  transport 
equation  to  exactly  recover  the  magnetogasdynamics  equations  of  (2-5).  Basically  there  is  no 
corresponding  “particle”  picture  to  represent  the  magnetic  field  evolution.  Therefore,  we  have  to  treat  the 
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flow  and  magnetic  fields  differently.  Whereas  the  flow  field  is  treated  with  the  gas-kinetic  theory,  the 
magnetic  fluxes  are  split  directly  based  on  the  macroscopic  equation  using  the  gas-kinetic  theory  as 
shown  in  [6]. 


The  governing  equations  of  (2-5)  can  be  rewritten  as 


8Q  dF  dG  dH  _ 
dt  dx  dy  dz 


(2-7) 


with 


Q= 


p 

pU 

pU 

pU2  +p~Txx 

pV 

pUV  -  Txy 

pW 

pUW-T „ 

< 

•,  F  =  ■ 

dT  r 

pe 

(pe  +  p)U  -  (i Ur „  +  Vr „  +  fVrX2)-  tc— 

Bx 

0 

By 

UBy-VBx 

Bz, 

UB2  -  WBr 

z  x 

(2-8) 


A  finite-volume  formulation  of  (2-7)  with  implicit  treatment  of  the  source  term  can  be  written  as 

a S_ 

8Q  J,  At  /=/,/,* 


[— - 
At 


f\n  _  r\n  77 «  in  _  t?  n  in 

1  jQ  _  *£/  &/  y  rl+\/2/il+V2  rl-\/2/il-\/2  j  ^n 

1  *  A(  L  T,  / 


(2-9) 


Here  F  is  the  numerical  flux  across  a  cell  interface  in  the  normal  direction 

pU 

pUU  +  pnx  -t„ 
pUV  +  pny  -  rny 
pUW  +  pnz 


F  = 


{pe  +  p)U  -  (1 Jr „  +  Vt  +  Wr „)  -  k 


ST 

ds„ 


(Vny+Wn2)Bx-U(Byny+B2n2) 


(Unx+Wnz)By-V(Bxnx  +Bznz ) 
(Unx+Vny)Bz-W(Bxnx+Byny) 


(2-10) 


where  U  =  Unx  +  Vn  +  Wnz ,  and  n  represents  the  normal  direction. 


Now  the  problem  that  remains  is  how  to  compute  the  numerical  flux  across  a  cell  interface.  As  mentioned 
earlier,  the  first  five  components  of  the  flux  of  (2-10)  are  related  to  the  flow  field.  Therefore,  they  can  be 
computed  by  the  gas-kinetic  BGK-NS  scheme  in  [7],  The  BGK  model  can  be  written  as 

ft+ufx+yf  +wfz=iG±  (2-11) 

T 

where  /  is  the  gas  distribution  function  and  g  is  the  equilibrium  state  approached  by  /  .  Both  /  and 
g  are  functions  of  space  (x,y,z),  time  t ,  particle  velocities  (u,v,w)  and  internal  degrees  of  freedom  g . 
The  particle  collision  time  r  is  related  to  the  viscosity  and  heat  conduction  coefficient.  The  equilibrium 
state  is  a  Maxwellian  distribution 
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f 


g=p\ 


m 


3+N 
1  2 


2  knT) 


exp{- 


m 

2kT 


[(u-U)2  +  (v-V)2  +(w-W)2  +  £2]} 


(2-12) 


where  m  is  the  molecular  weight,  k  is  the  Boltzmann  constant,  N  is  the  total  number  of  internal  degrees 

N 

of  freedom,  and  £2  =  Xg2  .  The  underlying  assumption  in  the  above  equilibrium  state  is  that  each  degree 

«=I 

of  freedom  shares  the  same  amount  of  internal  energy  kT  /  2  ,  or  the  so-called  equilibrium  flow. 

The  relation  between  the  first  five  components  of  the  conservative  variable  in  (2-8)  with  the  gas 
distribution  function  is 


where  y/  is  the  vector  of  moments 


qf  = 


v  = 


p 

pU 

pv 

pW 
ps  j 


=  \y/fdE 


1 

u 

V 

w 


—(u2  +v2  +w2  +£2) 

12 


(2-13) 


(2-14) 


and  dE  =  dudvdwd^  is  the  volume  element  in  the  phase  space  with  d<^  =  d^  ••■d%N .  Since  mass, 
momentum,  and  energy  are  conserved  during  the  particle  collisions,  /  and  g  satisfy  the  conservation 
constraint 


j(g~f)^dE  -0 


(2-15) 


at  any  point  in  space  and  time.  The  relation  between  the  first  five  components  of  the  fluxes  in  (2-8)  with 
the  gas  distribution  function  is  very  similar.  For  example, 


Ff 


pU 

pU2  +p-rxx 
pUV-r v 
pUW-T a 


f  =  J  y/ufdE 


(ps  +  p)U  ~(Utxx  +  Vt^  +  Wt^  )  -  k 


dT_ 
dx . 


(2-16) 


The  derivation  of  the  Navier-Stokes  equations  from  the  BGK  model  can  be  found  in  [8]. 


Using  (2-16),  the  first  five  components  of  the  fluxes  at  the  interface  can  be  computed  from  the  gas 
distribution  function  at  the  interface.  The  general  solution  of  the  BGK  model  gives  the  gas  distribution 
function  at  a  cell  interface  ( xM ,  2 ,  y j ,  zk  )  and  time  t  as 


f(xi+m,yJ,zk,t,u,v,w,4)  =  -lg(x',y',z,,t',u,v,w,4)e  t)/Tdt'  +  e  */„ (xM/2  - ut,y , -vt,zk  - wt ) 

To 
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(2-17) 

where  x'  =  xM/2  -  u(t  - 1') ,  y  =  y}  -  v(t  - 1')  ,  and  z'  =  zk  -w[t- 1’)  are  the  particle  trajectory  and  /0 
is  the  initial  gas  distribution  function  at  the  beginning  of  each  time  step  (t  =  0).  Two  unknowns  g  and 
f0  in  (2-17)  must  be  specified  in  order  to  obtain  the  solution  /  .  For  simplicity,  the  directional  splitting 
approach  is  adopted.  Therefore,  the  one-dimensional  notions  and  the  notion  of  (  *,+1/2  =  0>Ty  =0,Zk=0) 
will  be  used  in  the  following  text. 

Based  on  the  Chapman-Enskog  expansion  of  the  BGK  model,  the  gas  distribution  function  up  to  the 
Navier-Stokes  order  at  the  point  x  =  0  and  time  t  =  0  has  the  form 

ot  ox 


where  =  -r(—  +  u satisfies  the  compatibility  condition  \y/<f>xdE  =  0.  To  the  second-order 
accuracy,  the  gas  distribution  function  around  the  point  x  =  0  at  time  t  =  0  can  be  approximated  as 

ox  ot  ox 

Therefore,  given  the  initial  discontinuous  macroscopic  variables  at  the  left  and  right  hand  sides  of  a  cell 
interface,  the  initial  gas  distribution  function  f0  has  the  form 


\g,[\  +  a,x-r{alu  +  A,)\ 
°  \gr\\  +  arx-z{aru  +  Ar)] 


x<0 
x>  0 


(2-20) 


It  is  noteworthy  that  the  terms  proportional  to  r  in  (2-20)  represent  the  non-equilibrium  parts  of  the 
Chapman-Enskog  expansion.  These  non-equilibrium  terms  have  no  direct  contribution  to  the  macroscopic 
conservative  variables,  i.e., 


\{alu  +  Al)y/gldE  =  0 
\{aru  +  Ar)y/grd£  =  0 


which  are  also  the  equations  used  to  determine  A1  and  A r .  On  the  other  hand,  they  do  affect  the  fluxes. 
In  other  words,  the  gas-kinetic  BGK  scheme  has  more  information  and  gives  a  more  realistic  description 
of  the  viscous  flows  than  the  traditional  continuum  Navier-Stokes  solvers. 


Similarly,  one  can  construct  the  equilibrium  state  g  around  (x  =  0,t  =  0)  as 


g  =  &o[l +  0  -  H(x))a  x  +  H(x)a  x  +  At\ 


(2-22) 


where  g0  is  the  value  of  the  Maxwellian  distribution  function  at  ( x  =  0,  Z  =  0 )  and  H(x)  is  the  Heaviside 
function  defined  as 


H[x]  =  - 


1 


x<0 
x>  0 


(2-23) 


As  shown  in  the  following,  all  unknowns  of  a1 ,  aT ,  A1 ,  Ar ,  a1 ,  ar ,  and  A  in  both  (2-20)  and  (2-22) 
can  be  related  to  the  derivatives  of  a  Maxwellian  in  space  and  time. 


The  dependence  of  a1 ,  ar ,  A1,  Ar ,  a1 ,  a r ,  and  A  on  the  particle  velocities  can  be  obtained  from  a 
Taylor  series  expansion  of  a  Maxwellian  and  has  the  following  form 
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(2-24) 


a1  =a[  +  a'2u  +  a[v+  a[w+  a[  ~(w2  +  v2  +  w2  +£2)  =  ci,ay/a 
A1  =  A[  +  Al2u  +  A3v  +  A4w  +  A\  ~(u2  +  y2  +  w2  +  ^2)~  A,ay/a 


A  =  Ax  +  A2u  +  A3v+  A4w+  As  -^(u2  +  v2  +  w2  +  £ 2)-Aay/a 


where  all  coefficients  a[ ,  a'2 ,  Als  are  local  constants. 


After  the  reconstruction  stage,  one  obtains  the  left  and  right  hand  side  values  of  the  conservative  variables 
at  the  interface,  Qf+(,2  and  Q^n  ■  With  the  relation  of  (2-13),  the  values  of  p' ,  U1 ,  V! ,  W1 ,  T1  and 

pr ,  Ur ,  V ,  Wr ,  Tr  in  the  Maxwellian  distribution  functions  can  be  determined.  Similarly,  a1  and  ar 
can  be  computed  by 


ig'a'ysd--8^2  Q‘ 


*MI2  ~Xi 


\grary/d~=9M  Qi+'V2 


(2-25) 


xM  - x ; 


7+1/  2 


After  determination  of  a1  and  ar ,  A1  and  Ar  can  be  obtained  by  the  compatibility  condition  of  (2-21). 
Moreover,  according  to  [7],  the  values  of  p0 ,  U0 ,  V0,  W0 ,  and  ro  in  g0  can  be  determined  by 

\g0ydE  =  Ql  =\u>0\g'y/dE  +  \u<0fgry/dE  (2-26) 

Accordingly  a 1  and  ar  can  be  computed  by 


/g0aV-S’  =  — — — 
■*<+1/2  ~  X/ 


2°F 


(2-27) 


X; 


1+1  Xi  +  l/2 


The  only  unknown  left  is  A  ,  which  can  be  found  by  the  integration  of  the  conservation  constraint  of  (2- 
15)  over  the  whole  time  step  At 


\o‘l(S-f)y/dtdE  =  0 


(2-28) 


After  determination  of  the  distribution  function  /  at  the  cell  interface,  the  first  five  components  of  the 
fluxes  in  (2-8)  can  be  obtained  through  (2-16).  We  can  also  evaluate  the  heat  flux  across  the  cell  interface 

q  =  \\(u-  £/)[(«  -  U)2  +  (v  -  V)2  +  (w  -  Wf  +  £2  ]fds  (2-29) 

In  order  to  fix  the  problem  of  the  fixed  unit  Prandtl  number  for  the  BGK  solutions,  the  energy  flux  is 
modified  to  have  the  correct  heat  transfer  through  the  cell  interface,  F*/w  =  F p€  +  (l/Pr-l)g.  This 
completes  our  discussion  of  the  discretization  of  the  fluxes  related  to  the  flow  field. 
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Next,  we  discuss  how  to  compute  the  components  of  the  fluxes  related  to  the  magnetic  field.  Let  us  use 
the  flux  component  in  the  x-direction  as  an  example.  As  shown  above,  in  the  gas-kinetic  theory,  the  flux 
is  associated  with  the  particle  motion  across  a  cell  interface.  The  splitting  of  the  flux  component  in  the  x- 
direction  is  determined  by  the  particle  motion  in  this  direction.  Other  quantities  including  the  magnetic 
field  can  be  considered  as  passive  scalars  which  are  transported  with  the  x-direction  particle  velocity.  For 
example,  the  density  p  can  be  split  into 

p+  =  Lo \sdE = p\er^c^f^U) ~ p < u° >+  >p~ =  i<o isds = p\erfc^^u^ = P < u°  >-  (2-3°) 

Any  macroscopic  quantity  Z  without  explicitly  containing  the  x-component  velocity  U,  including  the 
magnetic  field  Bx ,  By,  and  Bz ,  can  be  split  in  a  similar  way 

Z+  =Z<m°  >+,  Z‘  =  Z<w°  >_  (2-31) 

The  above  relations  mean  that  the  quantity  Z  is  simply  advected  with  the  particle  transport  in  the  x- 
direction. 


On  the  other  hand,  the  x-direction  momentum  pU  can  be  split  into 


exp( — —U2) 

(PU)+  =l>0\ugd~=p[U<u°  >+  +- - =!M= - ]  =  p<ul  >+ 


m 

2 kT' 


-n 


(2-32) 


exp( — —  f/2) 

(puy  =  \u«*\u8ds  =  P\U  < u°  >-  -T - -ML= - ] ~p<u1  >. 

L  \m 

\2kTn 

Similarly,  any  quantity  containing  the  U  term,  including  UBX ,  JJBy ,  and  UBZ ,  can  be  split  as 

(zuy  =z<u]  >+ ,  ( zuy  =z<ux>_  (2-33) 

For  the  magnetic  field,  the  above  splitting  implies  that  the  field  is  frozen  into  the  particle  motion  and 
transported  with  the  fluid. 


As  a  result,  the  last  three  components  of  the  fluxes  in  (2-8),  which  are  related  to  the  magnetic  field,  can  be 
computed  by 


Fb  = 


UBy  -  VBX 
UB ,  -WB 


X  ) 


=(FBy  +(fb) 


B\- 


(2-34) 


The  positive  flux  is 


(FBy  = 


By  <«’  >+  —BXV  <  u°  >+ 
Bz  <ux  >+  -BXW <u°  >+  j 


(2-35) 


and  the  negative  flux  is 
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(2-36) 


t fbt  = 


0 

Bv  <m'  >_  -BXV  <u°  >_ 
B<u.'>  -BJW  <u°  > 


Finally,  we  will  discuss  how  to  numerically  enforce  the  divergence  free  condition  of  (2-6)  for  the 
magnetic  field.  A  common  approach  is  to  use  a  correction  method,  in  which  a  Poisson  equation  for  the 
potential  <f> 


d2d>}  dBj  A 
dx'dx1  8xJ 

is  solved  and  the  corrected  magnetic  field  B'c  is  obtained  through 


Bl  =B‘  + 


dx‘ 


where 


— r-  =  0  is  satisfied. 
dx' 


(2-37) 


(2-38) 
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3.0  ADJOINT  OPTIMIZATION  METHOD  FOR  MAGNETOGASDYNAMICS 
IN  THERMAL  AND  CHEMICAL  EQUILIBRIUM 


The  main  objective  of  this  program  is  to  develop  and  implement  viable  approaches  for  the  design 
optimization  of  magnetogasdynamics  problems.  To  further  clarify  this  objective,  we  intend  to  develop  a 
method  that  enables  us  to  predict  the  sensitivity/gradient  of  arbitrary  cost  (and  constraint)  functions  to 
very  large  numbers  of  design  parameters  with  reasonable  computational  cost.  At  the  moment,  we  have 
allowed  design  parameters  of  two  separate  kinds,  although  other  parameters  can  be  included  since  the 
formulation  is  entirely  generic: 

•  Vehicle  shape:  parametric  design  changes  to  a  baseline  geometry  that,  through  the  effect  of 
boundary  conditions  on  the  flow,  can  dramatically  change  the  functions  of  interest. 

•  Applied  electric  and  magnetic  fields:  the  strengths  and  shape  of  these  fields  can  have  a 
significant  impact  on  the  functions  of  interest.  Furthermore,  conductivity  and  permittivity 
distributions  can  also  impact  the  values  of  the  functions  of  interest  and  may  require  a  very  large 
number  of  parameters  to  be  appropriately  described. 

A  basic  requirement  of  our  approach  is  that  the  (cost)  functions  of  interest  (which  may  be  either  cost 
functions  for  the  optimization  or  constraint  functions  that  the  optimization  problem  must  satisfy)  can  be 
arbitrary  functions  of  the  flow  solution/state  variables.  In  our  early  work  we  have  tackled  the  drag 
coefficient  of  the  body  and  the  integrated  heat  flux  through  the  wall  as  basic  cost  functions,  but  the 
methodology  allows  for  any  other  function  to  be  included  with  a  very  small  development  cost,  namely, 
the  computation  of  the  partial  derivative  of  this  function  with  respect  to  the  state  vector  linearized  about 
the  solution  provided  by  the  flow  solver. 

With  this  in  mind  we  have  developed  a  method  that  can  easily  provide  derivatives  of  each  of  these 
functions  with  respect  to  arbitrary  numbers  of  design  parameters  (appropriately  chosen  by  the  designer) 
with  a  cost  which  is  equivalent  to  two  flow  solutions:  a  traditional  flow  solution  followed  by  an  adjoint 
solution  which  is  of  very  similar  cost  to  the  flow  solution  (in  terms  of  CPU  and  memory  requirements.) 
This  cost  is  for  the  derivative  of  a  single  cost  function  with  respect  to  a  large  number  of  design 
parameters.  If  the  gradients/sensitivities  of  additional  functions  are  required,  the  flow  solution  may  be 
reused,  but  additional  adjoint  solutions  (one  per  function)  must  be  carried  out  to  obtain  the  necessary 
gradients.  Note,  however,  that  the  treatment  of  additional  functions  does  not  require  the  development  of 
an  entirely  different  adjoint  solver  (this  misconception  is  common  in  our  community):  it  simply  requires 
that  a  new  right  hand  side  (one  could  call  this  a  forcing  vector)  be  developed.  As  mentioned  above  the 
cost  of  this  is  truly  minimal. 

After  careful  examination  of  the  governing  equations  of  the  problem  (see  Section  2.1)  the  choice  was 
made  to  develop  discrete  adjoint  equation  sets  for  the  N-S  low  magnetic  Reynolds  number 
approximation  in  Eq.  (26)  of  Ref.  [3].  This  choice  was  justified  due  to  the  extreme  similarity  of  the 
equations  with  the  Navier-Stokes  equations  of  fluid  flow,  for  which  we  have  built  up  considerable 
experience  through  the  years. 

Several  key  choices  in  the  development  have  also  been  made.  These  choices  make  this  program  unique  in 
comparison  with  previous  efforts  in  our  group  at  Stanford  University.  In  particular,  we  have  chosen  to: 

1 .  Use  the  discrete  adjoint  equations  (rather  than  the  continuous  formulation)  for  the  equations  of 
magnetogasdynamics. 
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2.  Develop  a  stand-alone,  unstructured,  discrete  adjoint  solver  for  the  magnetogasdynamics 
equations. 

3.  Choose  the  “discretization”  used  in  the  discrete  adjoint  formulation.  This  “discretization”  refers 
to  the  basic  discretization  for  which  the  adjoint  equations  have  been  developed. 

These  are  fundamental  choices  can  strongly  impact  the  usability  and  versatility  of  the  resulting  method 
and,  therefore,  details  on  the  rationale  for  these  choices  is  provided  below  after  a  short  description  of  the 
discrete  adjoint  method. 


3.1  Discrete  Adjoint  Equations 


The  conceptual  derivation  of  the  discrete  adjoint  equations  for  any  system  of  governing  equations  is 
straightforward  and  is  outlined  in  this  section.  We  consider  the  problem  of  the  minimization  of  an 
arbitrary  cost  function  of  interest,  J(U,(J),  which  involves  the  solution,  U,  to  the  a  particular  version  of 
the  magnetogasdynamics  equations  and  a  possibly  large  set  of  parameters,  /?,  that  may  include  various 
things  such  as  the  shape  of  the  body,  the  parametric  strength  of  imposed  magnetic  and  electric  fields,  and 
possibly,  the  location  and  strength  of  energy  addition.  The  changes  in  U  that  result  from  changes  in  /? 
are  not  arbitrary  since  they  must  satisfy  the  governing  equations  of  the  magnetogasdynamics  problem, 
N{U,/3)  =  N(U,X(J3))  =  0,  where  X  represents  the  mesh  that  may  (or  may  not)  change  as  a  result  of  the 
changes  in  the  design  parameters  /?  through  a  suitable  mesh  deformation  procedure.  The  total  derivatives 
of  the  cost  function  and  governing  equations  can  be  represented  as  follows 


If  we  solve  the  adjoint  equation 


dJ 

'  af 

dU 

dp_ 

.au. 

-dP- 

~cN 

_&J. 

dU 

U A 

+ 

-as  s> 

1 _ 1 

aj 

_|_  - 

vm 

=  0 


(3-1) 


dN 

T 

'a/  T 

YdUl 

.dU. 

,  where  v  is  the  adjoint  variable,  we  can  find  the 


derivatives  of  the  cost  function  as 


dJ 

oN 

VT  4- 

aj 

Up\ 

[ap\ 

V  -r 

SP- 

(3-2) 


Note  that  the  adjoint  equation  does  not  depend  on  [3.  It  only  has  to  be  solved  once  for  each  cost  function 


J .  The  derivatives  of  the  cost  function  in  Eq.  (3-2)  no  longer  involve  and  therefore  can  be 

calculated  without  re-computing  the  solution  to  the  governing  equations.  This  is  the  key  advantage  of  the 
adjoint  method.  For  systems  with  few  cost  functions  and  many  parameters/design  variables  the  discrete 
adjoint  approach  outlined  here  is  the  most  efficient  alternative  to  obtain  the  sensitivities  of  a  given 
function  of  interest.  Notice  that  all  of  the  expressions  above  are  discrete  versions  involving  vectors  of 
very  large  dimensions  and  very  large  and  sparse  matrices  (since  the  discretization  stencils  for  the  flow 
solver  tend  to  be  somewhat  compact.)  In  order  to  develop  the  discrete  adjoint  equations  one  needs  to 


come  up  with  expressions  for 


,dU\ 


and 


aj_ 


for  the  discrete  versions  of  N,  J ,  and  U. 


These 


expressions  must  take  into  account  the  actual  discretization  of  the  governing  equations  and  the  cost 
function.  Although  the  development  of  these  expressions  can  be  lengthy  and  laborious,  it  amounts  to 
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careful  application  of  the  chain  rule  of  differentiation  to  the  actual  discretization  stencil  chosen  for  the 
system  solution. 


All  of  the  expressions  for  the  adjoint  equation  (the  main  matrix, 


<5N 

-dU. 


with  the  boundary  conditions  that 


are  included  in  it  and  the  right  hand  side  vectors 


dJ_ 


for  the  two  cost  functions  described  earlier  have 


been  developed,  implemented  and  tested  (see  Section  3.4)  including  the  inviscid  and  viscous  terms  of  the 
residual  N,  a  first-order  discretization  of  the  artificial  dissipation  fluxes,  and  the  source  terms  for  the 
momentum  and  energy  equations  that  result  from  the  Lorentz  force  and  Joule  heating.  At  present,  the 
only  limitation  of  this  discrete  adjoint  implementation  is  that  the  variations  of  the  eddy  viscosity  with 
respect  to  the  flow  solution  have  been  neglected.  This  is  a  typical  approach  followed  when  the 
continuous  adjoint  formulation  is  used  as  it  is  impossible  to  analytically  derive  these  expressions.  In  our 
case,  it  is  simply  a  matter  of  additional  derivations  that  have  not  been  completed  yet. 


The  solution  of  the  adjoint  equation 


~sn 

T 

~  aj~ 

V  — 

-cU_ 

(3-3) 


is  achieved  using  a  matrix-free,  ILU(O)  preconditioned  GMRES  algorithm.  Since  Eq.  (3-3)  is  a  large 
linear  system,  GMRES  has  proven  to  be  one  of  the  most  robust  solvers  for  this  type  of  problem.  We  have 
chosen  to  use  a  matrix-free  version  (in  other  words,  we  never  store,  even  in  sparse  format,  the  matrix 


cN 

,dU. 


or  its  transpose  in  order  to  avoid  memory  storage  penalties.  Although  for  the  small  test  problems 


that  we  have  used  this  memory  penalty  is  not  truly  significant  and  does  not  hamper  our  efforts,  for  larger 
computations  it  may  become  a  significant  problem.  Instead,  we  recomputed,  on-the-fly,  the  terms  in 


„cU_ 


for  every  iteration  of  the  GMRES  algorithm. 


Alternative  compromises  between  storage  and  re¬ 


computation  (for  some  of  the  most  expensive  terms)  will  be  explored  in  Phase  II.  In  addition,  we  will 
enhance  the  performance  of  the  discrete  adjoint  solver  by  using  improved  preconditioning  techniques  and 


optimized  implementations  of  the  functions  that  calculate  the  various  Jacobian  terms  of 


(N_ 

,dU. 


3.2  Discrete  vs.  Continuous  Adjoint  Equations 

Since  most  of  our  efforts  at  Stanford  University  (Alonso,  Jameson,  Reuther,  et  al.)  have  focused  on  the 
development  of  continuous  adjoint  formulations  of  the  Reynolds-Averaged  Navier-Stokes  equations  it  is 
appropriate  to  discuss  in  this  section  the  rationale  that  led  us  to  the  choice  of  a  discrete  adjoint  approach 
for  the  equations  of  magnetogasdynamics.  The  various  choices  were  based  on  the  suitability  (and 
generality)  for  long  term  efforts  in  magnetogasdynamics. 

Several  key  reasons  have  led  us  to  the  decision  to  use  a  discrete  formulation  of  the  adjoint  equations: 


•  A  discrete  formulation  allows  for  the  choice  of  arbitrary  cost  functions  that  will  be  needed  to 
properly  pose  generic  magnetogasdynamics  optimization  problems.  In  contrast,  the  continuous 
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adjoint  approach  will  not  allow  proper  treatment  of  all  choices  of  cost  functions,  thus  limiting  the 
usability  of  the  resulting  method. 

•  The  derivation  of  the  continuous  adjoint  approach  that  includes  all  of  the  necessary  variations 
(turbulence  models,  electrical  conductivity,  permeability,  etc)  was  attempted  but  led  to  a  very 
high-level  of  complexity  unless  some  of  these  variations  were  neglected.  In  closely  coupled 
problems  in  magnetogasdynamics,  we  felt  that  neglecting  some  interactions  between  the  flow 
field  and  the  magnetic/electric  fields  would  lead  to  erroneous  results.  The  discrete  approach 
accounts  for  these  variations. 

•  The  derivatives  produced  by  the  discrete  method  are  in  agreement  with  the  discrete  computations 
of  the  system  solution. 

3.3  Discretization  Choices  for  the  Discrete  Adjoint  Equations 

The  resulting  discrete  adjoint  equations  depend  on  the  discretization  approach  chosen  for  the  residuals  of 

the  governing  equations.  Previous  experience  in  the  continuous  adjoint  area  ([28])  has  shown  that  one 
may  choose  consistent  discretizations  of  the  governing  equations  and  the  adjoint  equations  separately 
with  minimal  errors  for  reasonably-sized  meshes.  For  this  reason,  our  discrete  adjoint  solver  is  based  on  a 
finite-volume,  second-order  accurate  discretization  of  the  governing  equations  of  the  flow.  Although  this 
discretization  is  different  from  the  BGK  approach  used  in  the  flow  solver,  the  sensitivity  information 
must  be  consistent  in  the  limit  of  small  mesh  sizes.  Since  the  eigenvalues  of  the  flow  linearization  and  of 
the  adjoint  equations  are  identical,  the  only  requirement  on  the  choice  of  discretization  is  that  it  be  stable 
so  that  the  adjoint  solver  can  be  made  to  converge  ([27]). 

With  these  choices  for  the  discrete  adjoint  solver,  we  believe  that  we  are  in  a  position  to  tackle  the 
difficult  problem  of  magnetogasdynamics  design  regardless  of  the  flow  solver  and  discretization  used  and 
for  any  choice  of  cost  function  and  design  parameterization  that  can  involve  more  than  simple  changes  to 
the  physical  shape  of  the  geometry.  This  will  not  be  the  case  if  alternative  approaches  to  the  formulation 
of  the  adjoint  system  (for  example,  the  continuous  version)  are  employed. 

3.4  Sample  Sensitivity  and  Convergence  Histories  of  the  Discrete  Adjoint  Implementation 

The  procedure  outlined  above  has  been  implemented  into  a  stand-alone  discrete  adjoint  solver  for  the  low 
magnetic  Reynolds  number  approximation  of  the  magnetogasdynamics  equations.  As  mentioned  above, 
the  discrete  adjoint  solver  includes  all  of  the  necessary  terms: 

•  Inviscid  fluxes. 

•  Viscous  fluxes  (with  “frozen”  eddy  viscosity.) 

•  First-order  artificial  dissipation  fluxes. 

•  Source  terms  (for  both  the  Lorenz  force  and  Joule  heating.) 

operator,  we  have  developed  right  hand  side 

it  cost  functions  ,  J,  namely  the  coefficient  of  drag  of 

the  body  and  the  integrated  heat  transfer  through  the  (constant  temperature)  wall.  Additional  cost 
functions  can  be  easily  derived  and  implemented  as  another  right-hand-side  for  the  adjoint  solver. 


In  addition  to  these  terms  that  are  part  of  the 

T 


LtfC/J 


approximations  for  the 


aj 


IcU. 


terms  for  two  differe 
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The  ultimate  use  of  the  adjoint  solution  v  is  to  include  it  in  Eq.  (3-2)  to  compute  the  total  derivative  of  a 
cost  function  J  with  respect  to  the  design  parameters  /?.  The  ability  to  compute  these  derivatives 
accurately  is  the  main  validation  step  of  an  adjoint  solution  approach.  The  configuration  used  as  a 
preliminary  test  case  is  taken  from  Ref.  [3].  It  consists  of  a  blunt  cylinder  immersed  in  an  incoming  flow 
with  a  Mach  number  of  16  aligned  with  the  axis  of  symmetry  of  the  cylinder.  A  typical  flow  solution 
(obtained  with  a  finite  volume  implementation  of  the  low  magnetic  Reynolds  number  approximation  in 
the  parallel,  multiblock-structured,  Stanford  University  code,  TFL02000)  can  be  seen  in  Fig.  3.1  below. 


-i - 1 - 1 - 1 - ¥ 

X 


Figure  3-1  Mach  Number  Contours  for  the  Flow  Around  a  Blunt  Cylinder  at  a  Zero  Angle  of 
Attack.  Baseline  Solution  for  Discrete  Adjoint  Calculation. 

Based  on  the  baseline  solution  in  Fig  3-1,  we  computed  discrete  adjoint  solutions  for  both  the  coefficient 
of  drag  and  heat  transfer  cost  functions  described  earlier.  The  robustness  of  the  preconditioned  GMRES 
implementation  can  be  seen  in  Fig  3-2  below  where  the  convergence  history  of  the  discrete  adjoint 
solution  (with  J  being  the  coefficient  of  drag)  is  shown.  In  less  than  6000  residual  evaluations, 
corresponding  to  approximately  400  GMRES  iterations,  the  adjoint  equation  has  converged  over  3-4 
orders  of  magnitude  (the  various  graphs  denote  the  convergence  histories  of  the  adjoint  variables  for  the 
mass,  momentum,  energy,  and  turbulent  kinetic  energy  and  dissipation  rate  equations.  Previous  studies 
for  the  Navier-Stokes  equations  ([29],[30])  have  shown  that  this  convergence  level  in  the  adjoint  solution 
is  more  than  adequate.  In  fact,  a  level  of  convergence  of  2  orders  of  magnitude  has  been  found  to  be 
sufficient  for  accurate  gradient  information  to  be  obtained.  The  first  portion  of  the  convergence  history 
(labeled  I)  shows  the  convergence  for  a  mesh  that  has  been  coarsened  (by  taking  every  other  point)  in 
each  coordinate  direction.  After  a  certain  level  of  convergence  is  obtained,  the  fine  mesh  (128  x  96  x  16 
cells  in  the  streamwise,  normal,  and  azimuthal  directions)  is  started  (labeled  II  in  the  Fig.)  and 
convergence  proceeds. 
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Figure  3-2  Discrete  Adjoint  Solver  Convergence  for  the  Baseline  Solution  in  Fig  3-1  in  Both  the 

Fine  (II)  and  First  Level  Coarse  (I)  meshes. 


Using  Eq.  3-2,  we  can  use  this  adjoint  solution  (and  an  equivalent  one  obtained  for  J  being  the  heat 
transfer  measure)  to  compute  and  validate  the  gradients  with  respect  to  design  parameters  /?  of  interest. 
Following  the  types  of  design  parameters  that  we  have  mentioned  earlier,  we  present,  in  Table  3-1  below, 
the  results  of  the  comparison  of  the  discrete-adjoint-based  gradient  with  those  obtained  by  finite- 
differences  of  the  original  flow  solver.  Two  types  of  parameters  are  considered:  those  related  to  the 
strength  and  location  of  the  magnetic  field  and  those  related  to  the  shape  of  the  body  itself.  For  the 
magnetic  field  design  variables,  we  consider  two  design  parameters:  the  strength  of  a  dipole  and  its 
location.  For  the  shape  functions,  we  distributed  10  Hicks-Henne  ([30])  “bump”  functions  evenly  spaced 
from  the  stagnation  point  location  to  the  end  of  the  spherical  portion  of  the  body.  The  results  for  the  first 
and  last  Hicks-Henne  bump  functions  are  reported  here.  All  results  are  within  1-5%  of  the  finite 
difference  values  of  the  sensitivities  (appropriately  non-dimensionalized.)  We  are  currently  looking  into 
the  cause  of  these  differences  but  they  are  entirely  expected  as  the  exact  equivalence  can  only  be  obtained 
when  the  finite  difference  gradients  are  exact  (depending  on  the  choice  of  step  size)  and  when  all  of  the 

....  ,  .  ,  .  ,  r<wT 

terms  in  the  governing  equations  are  taken  into  account  when  computing  the  -  operator.  As 


VdUl 


mentioned  earlier  we  are  currently  using  a  first-order  approximation  to  the  artificial  dissipation  fluxes  in 
the  discrete  adjoint  solver  and  the  eddy  viscosity  field  is  “frozen”. 
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Cost  Functions,  J 

cD 

Integrated  Heat  Transfer 

Adjoint 

Finite  Diff. 

Rel.  Error  % 

Adjoint 

Finite  Diff. 

Rel.  Error  % 

Bump  1 

3.3227 

3.3789 

1.66 

37.3885 

35.6723 

4.81 

Bump  10 

0.3422 

0.3455 

0.95 

4.7523 

4.6676 

1.81 

Dipole 

Strength 

10.7588 

11.2377 

4.26 

26.3412 

27.0301 

2.55 

Dipole 

Position 

2.1577 

2.1553 

1.11 

7.5962 

8.0345 

5.45 

Table  3-1  Validation  of  Sensitivities  Using  Discrete  Adjoint  Solver 


We  are  currently  conducting  more  extensive  sensitivity  validation  numerical  experiments  and  providing 
the  computed  sensitivities  to  the  SNOPT  non-linear  constrained  gradient-based  optimizer  to  demonstrate 
the  use  of  our  adjoint  methodology  to  simple  design  problems.  We  expect  to  present  these  results  at  the 
time  of  our  final  presentation  at  the  end  of  June  2005. 
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4.0  BGK-BURNETT  SCHEME  FOR  GAS  DYNAMICS 
IN  THERMAL  AND  CHEMICAL  EQUILIBRIUM 


The  BGK-NS  solver  developed  above  may  be  insufficient  for  hypersonic  applications.  This  is  mainly 
because  a  space  vehicle  may  operate  in  a  high-speed  and  low-density  environment,  in  which  the 
rarefaction  effects  are  not  negligible.  A  key  measure  of  the  rarefaction  effects  is  the  Knudsen  number  K„, 
the  ratio  of  the  mean  free  path  to  the  characteristic  length.  The  higher  the  Knudsen  number  is,  the  more 
important  the  rarefaction  effects  become.  During  their  atmospheric  reentry,  the  freestream  Knudsen 
numbers  of  both  HITEN  and  Magellan  spacecrafts  range  from  2  to  6.3  [9],  indicating  that  the  flowfield  is 
in  the  transitional  regime.  In  the  flowfield  around  Reaction  Control  System  (RCS),  which  will  be  used  on 
future  Reusable  Launch  Vehicles  (RLV)  and  planetary  probes,  the  Knudsen  number  varies  from  much 
less  than  one  to  greater  than  one  from  the  continuum  jet,  through  the  interaction  region  and  into  the 
rarefied  freestream  [10].  The  accurate  simulation  of  these  mixed  continuum  and  rarefied  flows  calls  for  a 
computational  algorithm,  which  can  cover  all  continuum,  transitional  and  rarefied  flow  regimes,  or  at 
least  can  cover  the  Knudsen  number  up  to  one.  Such  a  need  is  further  emphasized  by  the  difficulty  of 
simulating  low-density,  high-enthalpy  flows  around  space  vehicles  on  the  ground-based  facilities. 

To  extend  our  gas-kinetic  scheme  beyond  the  continuum  regime,  instead  of  truncating  the  Chapman- 
Enskog  expansion  of  the  BGK  model  at  the  Navier-Stokes  order  in  (2-18),  we  keep  the  Chapman-Enskog 
expansion  of  the  BGK  model  up  to  the  Burnett  order 

/Burnet,  =  g  ~  *Dg  +  ZD(rDg)  (4-1) 

d  d 

where  D  =  —  +  u  — .  Since  we  are  going  to  develop  a  scheme  for  the  local  time  evolution  around  a  cell 
dt  dx 

interface,  the  variation  of  the  collision  time  r  around  a  cell  interface  within  a  time  step  is  ignored. 
Therefore,  the  equation  of  (4-1)  can  be  rewritten  as 

f  Burnett  =  g  ~  ^Dg  +  T* D* g  (4-2) 

where  <j>2  =  r2D2g  satisfies  the  compatibility  condition  j y/<f>2dE  =  0 .  Note  that  the  spatial  and  temporal 
variation  of  r  at  different  numerical  cells  and  different  time  steps  are  still  accounted  in  the  current  BGK- 
Bumett  scheme  by  changing  r  from  cell  to  cell  according  to  the  viscosity  coefficient  and  local  pressure. 


To  third-order  accuracy,  the  gas  distribution  function  around  the  point  x  =  0  at  time  t-  0  can  be 
approximated  as 


fi 


Burnett 


dg  1  32g  2 

:g  +  —  x  + - f-jc 

dx  2  dx2 


■r <&  +  „&)■ 

dt  dx 


dtdx  dx2  dt 2  dtdx  dx2 


)  (4-3) 


Therefore,  given  the  initial  discontinuous  macroscopic  variables  at  the  left  and  right  hand  sides  of  a  cell 
interface,  the  initial  gas  distribution  function  /„  has  the  form 


gl{\  +  alx  +  ^b'x2  -T[a'u  + A'  +(C!  +b'u)x]  +  T2(B'  +  2uC!  +  b'u2)},x< 0 
gr{\  +  arx  +  ^brx2 -v[aru  +  Ar  +(Cr  +bru)x]  +  T2(Br  +  2uCr  +bru2)},x>  0 


(4-4) 


The  values  of  p1 ,  U1 ,  V1 ,  W' ,  T1  and  pr ,  Ur ,  Vr ,  Wr ,  Tr  in  the  Maxwellian  distribution  functions, 
and  a1 ,  ar  are  computed  in  the  same  way  as  in  the  BGK-NS  scheme.  The  new  unknowns  bl  and  br  can 
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be  determined  from 


\  g'b'ydS  =  &-)', M , 

&2  (4-5) 

j  g'b'Vds=<?£yMn 

OX 


After  determination  of  a1 ,  ar  and  b! ,  br ,  A1  and  Ar  are  still  obtained  by  the  compatibility  condition  of 
(2-21).  The  new  imknowns  C1  and  Cr  can  be  determined  by  the  compatibility  condition 

\(b'u  +  C')y/gld~=  0 

(4-6) 

\(bru  +  Cr)ipgrdS  =  0 

and  B1 ,  Br  further  by  the  compatibility  condition 

l(b'u2+2uC'  +Bl)Wg‘dS  =  0 

\(bru2 +2uCr +Br)y/grdS  =  0  1  ' 


Similarly,  one  can  construct  the  equilibrium  state  g  around  (x  =  0,/  =  0)as 

g  =  go[i  +  (l  -  H(x))(a'x  +  i*V)  +  H(x)(7x  + 1 brx2)  +  Ai  +  ~Bt 2  +  Cxt ]  (4-8) 

The  unknowns  in  (4-8)  can  be  determined  in  a  similar  way  as  above. 

Next,  we  first  test  our  BGK-Bumett  scheme  for  the  Poiseuille  flow  under  the  external  forcing  term  with 
the  Knudsen  number  Kn  =0.1.  As  pointed  out  by  many  authors  [11-14],  even  for  this  simple  case  with 
relative  small  gradient  and  Knudsen  number,  the  Navier-Stokes  equations  fail  to  predict  qualitatively 
correct  solution.  Specifically,  for  the  external  force  driven  case,  the  Navier-Stokes  equations  fail  to 
reproduce  the  central  minimum  in  the  temperature  profile  and  a  non-constant  pressure  profile  in  the  cross¬ 
stream  direction,  which  are  both  predicted  by  the  gas  kinetic  theory  and  observed  in  the  DSMC  results. 
Furthermore,  based  on  the  Navier-Stokes  equations,  it  is  impossible  to  correct  this  failure  by  modifying 
the  equation  of  state,  transport  coefficient  or  boundary  conditions.  Unlike  the  slip  phenomena,  the 
discrepancy  is  not  just  near  a  boundary  but  throughout  the  system.  The  similar  discrepancy  is  also 
happening  for  the  pressure-driven  case. 


The  setup  of  the  external  force  driven  case  is  given  as  follows  [15-16].  The  simulation  fluid  is  a  hard 
sphere  gas  with  particle  mass  m- 1  and  diameter  d  = 1 .  At  the  reference  density  of  p0=  1 .21  x  10~3 ,  the 

mean  free  path  is  l0  =-  m 


-J2np0d 

temperature  is  T0  =  1 .  The  reference  fluid  speed  is  U0 


=  186.  The  distance  between  the  thermal  walls  is  Ly  =10/0  and  their 


2kTr 


m 


—  =  1 ,  so  the  Boltzmann  constant  is  taken 


as  k  -  — .  The  reference  sound  speed  is  c0  = 


=0.91  with  /=—  for  a  monatomic  gas.  The 
m  3 


p  kT 

reference  pressure  is  p0  — -  =  6.05  x  10  4 .  The  acceleration  and  pressure  gradient  are  chosen  so  that 


m 


the  flow  will  be  subsonic,  laminar,  and  of  similar  magnitude  in  the  two  cases.  Specifically, 
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A>|/|  =  8.31 x  10  8  for  the  acceleration-driven  case  and  =  1 .08  x  10  7  for  the  pressure-driven  case 

( p+  -  \.5p0,  p_  =  0.5 p0,  Lx  -  30/0 ).  In  both  cases,  the  Rnudsen  number  is  0.1  and  the  Reynolds  number 
is  of  order  one.  In  all  calculations,  the  cell  size  is  half  the  mean  free  path  of  the  initial  data. 


Figure  4- 1(a)  presents  the  results  for  the  force-driven  case  in  the  cross-stream  direction  using  the  above 
BGK-NS  scheme  with  the  slip  kinetic  boundary  condition.  The  circles  in  Figure  4-l(a)  are  the  well- 
verified  DSMC  results  [16].  Although  the  BGK-NS  scheme  is  an  accurate  Navier-Stokes  solver,  even 
with  the  slip  boundary  condition,  the  predicted  pressure  distribution  is  constant  in  the  cross-stream 
direction,  which  is  different  from  the  DSMC  solution.  In  [16],  a  different  Navier-Stokes  solver  is  used. 
Both  Navier-Stokes  solvers  give  qualitatively  similar  results.  In  order  to  resolve  the  discrepancy  between 
the  Navier-Stokes  and  DSMC  solutions,  the  gas-kinetic  BGK-Bumett  scheme  is  used.  Figure  4- 1(b) 
presents  the  BGK-Bumett  results.  It  is  found  that  the  curved  pressure  distribution  in  the  cross-stream 
direction  is  captured  well.  So,  up  to  the  Burnett  order,  the  non-constant  pressure  distribution  can  be 
obtained.  This  is  consistent  with  the  analysis  in  [13]. 


(a)  BGK-NS  vs  DSMC 


(b)  BGK-Bumett  vs  DSMC 

Figure  4-1  Flow  distribution  in  the  cross-stream  direction  for  the  external  force-driven  case 

(o  DSMC  results,  —  BGK  results) 
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In  order  to  further  validate  the  BGK-Bumett  approach,  we  have  calculated  the  mass  flow  rate  through  the 
channel  for  the  pressure-driven  Poiseuille  flow.  Both  BGK-NS  and  BGK-Bumett  schemes  are  used  in  the 


current  calculations.  The  normalized  mass  flow  rates  are  defined  by 


g 

pU'h 


,  where  g  is  the  mass  flow 


rate  and  pU  h  is  the  normalization  factor.  The  velocity  U  is  defined  by  U  =  a  ,  where  a  is 


related  to  the  pressure  gradient  in  the  channel  p  =  p0  (1  +  — -)  and  h  is  the  channel  width.  The  calculated 

h 

solutions,  as  well  as  the  analytical  solution  of  [17],  are  shown  in  Figure  4-2.  From  this  figure,  we  can 
clearly  see  the  improvement  of  the  BGK-Bumett  solution  over  the  BGK-NS  solution.  For  most  micro- 
channel  flows  in  the  laboratory,  such  as  the  experiments  in  [18],  the  highest  Knudsen  number  at  the  outlet 
is  about  0.2.  Therefore,  in  the  flow  regime  of  these  experiments,  the  slip  Navier-Stokes  equations  are 
capable  of  capturing  the  physical  solution.  However,  as  the  Knudsen  number  increases,  such  as  up  to  0.5, 
the  BGK-Bumett  scheme  should  be  a  more  appropriate  numerical  method  than  the  BGK-NS  solver. 


io  m 1  nr  in' 

Kn 


Figure  4-2  Relation  of  the  normalized  mass  flow  rate  versus  Knudsen  number  for  Poiseuille  flow 
(o  Boltzmann  solution,  x—  BGK-NS,  +  BGK-Burnett) 
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5.0  GENERALIZATION  OF  GAS-KINETIC  NS  SOLVER  TO  HYPERSONIC  FLOWS 

IN  THERMAL  AND  CHEMICAL  NON-EQUILIBRIUM 


Next,  we  discuss  how  to  generalize  the  above  BGK-NS  scheme  to  hypersonic  flows  in  thermal  and 
chemical  non-equilibrium.  As  a  first  step,  a  simpler  Kinetic  Flux  Vector  Splitting  (KFVS)  scheme  is 
constructed  here  for  the  governing  equations  adopted  by  a  widely  used  hypersonic  CFD  code,  LAURA. 


According  to  [1-2],  the  governing  equations  solved  by  LA  URA  can  be  written  as 

_d_ 
dxj 


^(pUi)  +  ^7(pUiUj  +  Sij p  -  TiJ)  =  0 
ot  dxJ 

!•(/*)  +  7T [<j*  +  P)Uj  -  pi(hsDs  ^-)  - U'T*  ~k~~~kv  ~~t)  =  0 
ot  dxJ  s= i  dxJ  8xJ  dx> 

i  (psv  )  +  TT  [P£vuJ  -  P  I  &r*Dt^)-Ky^]  =  &y 

ot  dxJ  s=mol.  dxJ  dxJ 


(5-1) 


Here  ps,  Ds ,  ys,  d>s ,  hs ,  and  hv  s  are  the  density,  effective  diffusion  coefficient,  mole  fraction,  the 
mass  production  rate  due  to  chemical  reactions,  enthalpy  per  unit  mass,  and  vibrational-electronic 
enthalpy  per  unit  mass  for  species  s  respectively.  A  1 1 -species  air  model  (N,  O,  N2,  02,  NO,  N+,  O4,  N2 , 

02 ,  NO+,  and  e  )  is  used  in  LAURA.  A  common  velocity  U‘  is  assumed  for  all  species.  To  account  for 
thermal  non-equilibrium,  a  two-temperature  model  is  used,  in  which  the  translational  and  rotational 
energy  modes  are  assumed  in  equilibrium  at  the  translational  temperature  T ,  and  the  vibrational, 
electronic,  and  electron  translational  energy  modes  are  assumed  in  equilibrium  at  the  vibrational 

temperature  Tv .  p,  p ,  r,J ,  s ,  ev  ,  k  and  kv  are  the  density,  pressure,  shear  stress,  total  energy  per 
unit  mass,  vibrational-electronic  energy  per  unit  mass,  frozen  thermal  conductivity  for  translational- 
rotational  energy,  and  frozen  thermal  conductivity  for  vibrational-electronic  energy  for  the  gas  mixture 
respectively.  cov  represents  the  vibrational-electronic  energy  production  rate  due  to  particle  collisions.  A 

detailed  description  of  the  governing  equations  of  (5-1),  including  thermodynamic  relations  and  chemical 
kinetic  models,  is  referred  to  [1-2], 


A  finite-volume  formulation  of  (5-1)  with  implicit  treatment  of  the  source  term  is  the  same  as  (2-9)  in  the 
form.  The  only  problem  left  is  how  to  construct  KFVS  scheme  for  computation  of  the  numerical  fluxes  of 
(5-1)  at  the  interfaces.  In  order  to  recover  the  governing  equations  of  (5-1),  1 1  BGK  models  are  needed 
such  that  each  species  satisfies  its  own  BGK  model 


dfs  t  dfs  =  Ss  ~fs 

ot  S  dxj  Ts 


For  flows  in  thermal  equilibrium,  the  Maxwellian  distribution  function  can  be  written  as 


(5-2) 


3 +N, 


Ss  =Ps 


m. 


2knT, 


S  J 


Wl  3  .  _  o 

exp{-— ML« -U\f  +I(<f/)2]} 

ZKl  j  j 


(5-3) 


which  is  basically  the  application  of  (2-12)  to  each  species.  For  hypersonic  flows  in  thermal  non¬ 
equilibrium  with  the  two-temperature  model  used  in  [1],  however,  one  has  to  generalize  the  Maxwellian 
distribution  function  to 
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&  s  Ps 


2k,7rTc 


exP{-^=-[S(«/  ~U{)2  +  Z(£/)2] 
2kTs  j  j 

m  N*  ■  w-£  ■  , 


2*7V,s  /  J 

where  £,{ ,  rjj ,  are  rotational,  vibrational,  and  electronic  degrees  of  freedom  respectively,  and  N* , 
Nvs  ,  Nf  are  the  numbers  of  rotational,  vibrational,  and  electronic  degrees  of  freedom  respectively.  For 
electrons,  Ts  =  TV  s . 

The  relation  between  the  macroscopic  variables  like  the  mass  ps ,  the  momentum  pUl ,  the  total  energy 
ps ,  the  vibrational-electronic  energy  psv  and  the  gas  distribution  function  fs  is 


pU‘  n 

Q=  =1  WJ'dS, 

ps  j=l 


where  the  moments  of  y/s  are  given  as  follows 


•  3  (UJ)2  N*  <pj\7  w"  (nj)2  Nf(pj\2  N';  (rjj)2  n?  (pj)2  „ 

Vx  =  (i,o,o,o,o,o,o,o,o,o,o  ,«;,E^r+  +  1^-+  l^r-,£^r-+  Z^r-)r 

j= i  2  j= l  2  7=1  2  y=i  2  7=i  2  7=1  2 

•  3  ('«-'■') 2  NUPJ\2  Kfrjj)2  N?(£j)2  Nl  (VY  Nf(pj)2  „ 

Vi  =  (0,1,0,0,0,0,0,0, 0,0,0, i^-  +  z^-+  z^-+  z^-)r  (5-6) 

y=i  2  y=l  2  7=1  2  7=1  2  7=1  2  7=1  2 


Due  to  the  momentum  and  energy  exchange  during  particle  collisions  between  all  species,  however,  the 
Maxwellian  distribution  functions  in  (5-4)  are  not  completely  independent.  Consistent  with  the  physical 
model  used  for  (5-1),  we  assume  gs  to  have  the  common  velocity  and  temperatures 


8  s  =  Ps 


3 +N’ 

m,  2  (  m. 


2knT  J  2 kjffr 


2  171  3  .  N?  .  rtf  Nvt  .  Nf  . 

exp{--^[S(«/  -C/O2  +  Z(^)2]-~^[Z(/7/)2  +  Z(C/)2]} 

2kT  j  j  2 kTy  j  j 


Using  (5-5),  one  can  find 


i  psui 

U '■  =  - - 

p 

2Zps(ss -Sy^-ipQJ1)2 -3pe—Ty 
r  =  ^ 

Z(i\Cj,+3)/7J*/mJ 


^=TT 


2 


Z(AA;  +N?)p,k/m, 
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The  relation  between  the  fluxes  of  (5-1)  and  the  gas  distribution  function  fs  is 


F  = 


psUJ-pD—ys 

dxJ 

TTtTTj  x  ,dU‘  dU\  2  3Uk  . 

W+pS,-*— 


TTTTi  dT  dTy  n  dy 
PHU  ~ri—-riv—--pY.hsD, 
dxJ  dxJ 


Tji  ,dU'1  .  dU\  .  2  TTi  dUk  e 

j'-'j  •  U  Z/( - : — \ - — )  ■+ - U  p - —  Su 

s  s  dxJ  dxJ  dx‘  3  dxk  J 


P£VUj  -Vyd-^-plhvsDs^ 
8xJ  tx  ,s  dxJ 


=  'LWMfsdZl 

s= 1 


(5-9) 


where  /  is  the  Chapman-Enskog  expansion  up  to  the  Navier-Stokes  order  like  the  one  given  in  (2-18). 


KFVS  scheme  can  be  simply  constructed  by  splitting  the  flux  term  of  (5-9)  into  the  positive  and  negative 
parts  according  to  the  sign  of  the  particle  velocity 


Fj+i/2  =  hl>oYs»Jsgtd~s -^~(L>0  V'suJsgfd^s)-T^J(lJ>0xps(uJs)2gkdSs) 
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+  lu!<oV/suigfdFs  -  rjt(\u,<0y/suJsg*dZs)-  r^T(^<o  VsW)2  gfdS,)] 


(5-10) 


The  time  derivative  term  in  (5-10)  can  be  replaced  by  the  spatial  derivative  terms  using  the  following 
compatibility  condition 
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The  thermodynamic  relations  and  chemical  kinetic  models  used  in  the  above  KFVS  scheme  follow  the 
work  of  Gnoffo  et  al.  in  [1-2]. 


To  test  this  non-equilibrium  KFVS  scheme,  we  simulate  the  RAM-C  II  flight  test  case.  The  configuration 
is  a  spherically  blunted  cone  with  the  nose  radius  of  0.1524  meter  and  the  half  cone  angle  of  9°.  The  total 
length  is  1 .3  meter.  The  ffeestream  speed  is  7650  m/s.  The  test  case  we  selected  is  the  one  at  the  altitude 
of  61  km.  As  a  result,  the  ffeestream  Mach  number  is  23.9  and  the  ffeestream  Reynolds  number  based  on 
the  nose  radius  is  about  1.95xl04. 


Figure  5-1  presents  the  computational  grid  used  for  the  computation,  which  covers  half  the  body  for  less 
computational  costs.  There  are  51  points  along  the  streamwise  direction,  61  points  along  the  azimuth 
direction,  and  61  points  along  the  normal  direction  with  the  first  grid  spacing  of  0.00026  nose  radii. 
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Figure  5-1  Computational  grid 


Figure  5-2  presents  the  predicted  nondimensionalized  pressure  {p / pm V!X2 )  contours  in  the  first  azimuth 
plane  with  the  fixed  wall  temperature  of  1500°K.  The  chemical  kinetic  models  used  are  from  [19-20]. 


The  predicted  translational-rotational  and  vibrational-electronic  temperatures  are  presented  in  Figure  5-3. 
The  results  are  consistent  with  the  results  obtained  in  [21-22],  The  highest  translational-rotational 
temperature  after  the  shock  around  the  nose  can  be  around  20,000°K  whereas  the  highest  vibrational- 
electronic  temperature  can  be  around  10,000°K.  In  such  a  high  temperature  environment,  the  flow  can  be 
weakly  ionized. 


0.15 


0.05 


0.15 


(a)  translational-rotational  temperature 


Figure  5-3  Translational-rotational  and  vibrational-electronic  temperatures 


10  11 

Figure  5-4  further  presents  the  total  ion  mass  fraction  (I  (ps  I  m5)l  t.(ps  I  ms))  contours  and  the  free 

s=6  5=1 

electron  mass  fraction  ((pe  / me)I  £(ps  / ms))  contours.  It  is  found  that  under  the  above  high 

5=1 

temperatures,  the  ionized  particles  account  for  only  about  0.7%  of  the  total  particles.  Such  a  weak 
ionization  is  apparently  not  enough  for  electromagnetic  flow  control.  More  importantly,  Figure  5-4 
indicates  that  the  total  ion  mass  fraction  is  almost  the  same  as  the  free  electron  mass  fraction.  Therefore, 
the  first  term  of  the  Lorentz  forces  added  to  the  momentum  equations  by  Appleton  and  Bray  in  [23], 

( nio„s  ~  ne)e(E‘  +  BkUj  -  BjUk ) ,  is  almost  zero.  Here  nions  is  the  total  number  density  of  the  ions,  ne 
is  the  number  density  of  free  electrons,  and  e  is  the  electron  charge.  The  electromagnetic  effects  can  be 
exhibited  only  through  the  conduction  current  density:  nee(U‘  -£/'). 
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(a)  total  ion  mass  fraction 


X 

(b)  electron  mass  fraction 

Figure  5-4  Ions  and  electron  mass  fractions 


In  the  current  literatures,  the  conduction  current  density  is  often  calculated  by  the  Ohm’s  law: 
j‘  -  <y(E '  +  BkUJ  -BJUk) .  According  to  the  definition  of  a  in  [24],  we  have  computed  the  electric 
conductivity  for  the  above  RAM-C  II  case  and  presented  the  predicted  contours  in  Figure  5-5.  The  results 
are  comparable  with  those  given  in  [24], 
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6.0  FUTURE  WORK  FOR  PHASE  n 


As  indicated  in  Figure  6-1,  which  is  taken  from  [25],  the  cylinder  and  sphere  can  be  considered  as  simple 
models  of  both  the  wing  leading  edge  and  the  nose  of  space  vehicles  respectively.  Therefore,  in  Phase  II, 
we  will  first  extend  the  BGK-Bumett  scheme  to  the  curvilinear  coordinates  and  use  the  simpler  cylinder 
case  to  examine  the  validity  Knudsen  number  range  of  the  BGK-Bumett  scheme.  The  test  data  in  [26] 
will  be  used  to  validate  the  BGK-Bumett  scheme. 


Figure  6-1  Similarity  between  simple  geometries  and  wing  leading  edge/nose  of  space  vehicles  [25] 

After  the  validation,  we  will  extend  the  BGK-Bumett  scheme  to  magnetogasdynamics  in  thermal  and 
chemical  non-equilibrium  and  develop  an  adjoint  optimization  method  based  on  this  solver  using  a 
similar  approach  to  the  one  outlined  earlier  in  this  document  and  ensuring  exact  discrete  consistency 
between  the  flow  and  adjoint  solvers  on  arbitrarily  shaped  geometries  and  flow  conditions. 
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